We  introduce  a  numerical  procedure  for  the  construction  of  interpolation  and  quadra¬ 
ture  formulae  on  bounded  convex  regions  in  the  plane.  The  construction  is  based  on 
the  behavior  of  spectra  of  certain  multiplication  operators  and  leads  to  nodes  which 
are  inside  a  prescribed  convex  region  in  The  resulting  interpolation  schemes  are 
numerically  stable  and  the  quadrature  formulae  have  positive  weights  and  almost  (but 
not  quite)  optimal  numbers  of  nodes.  The  performance  of  the  algorithm  is  illustrated 
by  several  numerical  examples. 
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1.  Introduction 


Numerical  integration  and  interpolation  constitute  one  of  oldest  areas  of  numerical 
analysis,  and  find  applications  in  most  branches  of  applied  mathematics.  In  general 
outline,  the  formulation  is  as  follows.  Given  a  region  G  G  (usually,  with  k  <  3) 
and  a  collection  of  functions  (pi,  ip2,  ■  ■  ■  ,^n  mapping  G  to  we  would  like  to  hnd  a 
collection  of  points  Xi,  0:2, . . . ,  such  that  for  each  point  x  in  G,  there  exist  coefficients 
Ui(t),  0:2 (t),  ... ,  Q:m{x)  such  that 

m 

(1.1)  f{x)  =  ^aj  f{xj), 

i=i 

for  all  /  :  G  — of  the  form 

n 

(1-2)  f{x)  =  ^j3iPi{x)- 

i=\ 

in  other  words,  the  interpolation  formula  (1.1)  is  exact  for  any  /  that  is  a  linear  com¬ 
bination  of  the  functions  {pt}. 

In  practice,  (1.1)  is  not  sufficient,  since  if  the  resulting  interpolation  formulae  are  to 
be  of  any  practical  use,  the  coefficients  {o;*}  must  not  be  very  large,  and  have  to  be  easy 
to  evaluate.  The  following  theorem  (see  [9])  guarantees  the  existence  of  numerically 
acceptable  interpolation  coefficients  under  extremely  weak  conditions. 

Theorem  1.1.  Suppose  that  S  is  an  arbitrary  set,  n  is  a  positive  integer,  /i,/2,  ■  ■  ■  ■,  fn 
are  bounded  complex-valued  functions  on  S,  and  e  is  a  positive  real  number  such  that 
e  <  1.  Then,  there  exist  n  points  Xi,X2,  ■  ■  ■  ,Xn  in  S  and  n  functions  gi,g2, . . .  ,gn  on  S 
such  that 

(1.3)  \9k{x)\  <  I e 

for  all  X  in  S  and  k  =  1,2, . . .  ,n,  and 

n 

(1-4)  f{x)  =  ^f{xk)gk{x) 

k=l 

for  all  X  in  S  and  any  function  f  defined  on  S  via  the  formula 

n 

(1-5)  f{x)  =  ^Cifi{x) , 

i=\ 

for  some  complex  numbers  Ci,  C2, . . . ,  c„. 

Thus,  the  existence  of  a  stable  interpolation  formula  is  never  a  serious  issue;  this 
paper  deals  with  numerical  procedures  for  constructing  such  interpolation  formulae. 


Numerical  integration  of  functions  is  closely  related  to  their  numerical  interpolation, 

but  involves  a  number  of  subtleties.  Given  a  region  G  G  and  a  collection  of  functions 
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Lpi,ip2)  ■  ■  ■  '■  ^  ^  one  wishes  to  find  points  Xi,X2,  ■  ■  ■  ,Xm  G  and  coefficients 

Wi,W2, . . . ,  Wm  e  such  that  for  each  of  the  functions  {v?i}, 


the  pair  {{a;*},  {wi}},  is  frequently  referred  to  as  a  quadrature  formula  for  the  functions 
{99*},  with  {xi}  the  nodes,  and  {tCj}  the  weights  of  the  quadrature. 

Again,  (1.6)  is  not  sufficient  for  the  quadrature  {{xj},  {tCi}}  to  be  of  practical  use. 
Indeed,  its  numerical  stability  will  be  poor  if  the  weights  {tCj}  are  large.  Obviously, 
if  all  of  the  weights  {wi}  are  positive,  this  issue  does  not  arise,  and  quadratures  with 
positive  weights  are  generally  preferred  (though  in  most  cases  not  absolutely  necessary). 

Obviously,  given  an  interpolation  formula  of  the  form  (1.4)  for  a  set  of  functions 
/i,  /2, . . . ,  /n,  it  is  easy  to  construct  a  quadrature  that  is  numerically  stable;  this  obser¬ 
vation  is  formalized  in  the  following  corollary  to  Theorem  1.1  (see  [9]). 

Corollary  1.1.  Suppose  that  S  is  a  measure  spaee,  w  is  a  nonnegative  real-valued 
integrable  funetion  on  S  (that  serves  as  the  weight  for  integration),  n  is  a  positive 
integer,  /i,/2,  ■  ■  ■  ,fn  are  bounded  complex-valued  integrable  functions  on  S,  and  e  <  1 
is  a  positive  real  number. 

Then,  there  exist  n  complex  numbers  Wi,W2,...,Wn  such  that 
(1.7)  M<(l+e)/u.(x)<ix 

for  all  k  =  1,2, ...  ,n,  and 

/n 

f{x)w{x)  dx  =  E  Wkf{Xk) 

k=i 

for  any  function  f  defined  on  S  via  the  formula  (1.5),  where  xi,X2, . . .  ,Xn  are  the  n 
points  in  S  chosen  in  Theorem  1.1. 

One  fundamental  difference  between  interpolation  and  integration  is  the  fact  that 
interpolation  in  an  n-dimensional  space  always  requires  at  least  n  points  (and,  according 
to  Theorem  1.1,  n  points  are  always  sufficient);  on  the  other  hand,  a  quadrature  formula 
for  integrating  n  functions  often  requires  fewer  than  n  points.  The  reader  is  referred  to 
Section  2.4  for  a  somewhat  more  detailed  discussion  of  this  class  of  issues. 

In  one  dimension,  the  existing  theory  of  quadrature  formulae  is  quite  satisfactory. 
Gaussian  quadratures  (see,  for  instance,  [3])  are  optimal  for  integrating  polynomials  up 
to  a  fixed  degree,  while  generalized  Gaussian  quadratures  are  used  for  integrating  more 
general  classes  of  functions  including  smooth  functions  and  functions  with  end-point 
singularities  [8,  12]. 

In  dimensions  higher  than  1,  much  less  is  known  regarding  optimal  numerical  integra¬ 
tion  formulae.  Classical  expositions  on  constructing  quadratures  for  multiple  integrals 
can  be  found  in  [4] ,  [6] ,  [7] ,  and  a  more  recent  overview  of  the  development  of  the  theory 
is  given  in  [2]. 

This  paper  is  partly  devoted  to  constructing  quadratures  on  convex  regions  in  two 

dimensions  exhibiting  Gaussian  behavior  (the  number  of  nodes  of  the  quadrature  is 
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smaller  than  the  number  of  functions  to  be  integrated)  using  a  simple  scheme,  which 
is  relatively  easy  to  implement.  The  quadratures  respect  the  symmetry  of  the  region, 
have  positive  weights  and  all  their  nodes  are  interior  to  the  region  of  integration. 

The  structure  of  this  paper  is  as  follows.  Section  2  introduces  numerical  preliminaries. 
In  Section  3  we  develop  the  mathematical  apparatus  to  be  used  by  the  algorithm  in 
Section  4.  We  summarize  the  numerical  results  in  Section  5. 


2.  Numerical  Preliminaries 

2.1.  Singular  Value  Decomposition.  The  singular  value  decomposition  (SVD)  is  a 
ubiquitous  tool  in  numerical  analysis.  In  the  case  of  real  matrices,  the  SVD  is  given  by 
the  following  lemma  (see,  for  instance,  [5]  for  more  details). 

Lemma  2.1.  For  any  nxm  real  matrix  A,  there  exist  a  positive  integer  p,  an  nxp  real 
matrix  U  with  orthonormal  columns,  anmxp  real  matrix  V  with  orthonormal  columns, 
and  a  pxp  real  diagonal  matrix  S  =  {sij)  whose  diagonal  entries  are  positive,  such  that 
A  =  U SV*  and  sa  >  for  all  i  =  1, . . .  ,p  —  1. 

The  diagonal  entries  su  of  S  are  the  nonzero  singular  values  of  A.  The  columns  of  U 
span  the  column  space  of  the  matrix  A  and  hence  provide  an  orthonormal  basis  of  the 
column  space  of  A.  This  property  of  the  SVD  is  used  in  Section  4  for  the  computation 
of  orthonormal  basis  functions. 


2.2.  The  QR  Method  for  fiudiug  the  spectrum  of  a  matrix.  There  are  several 
methods  available  for  determining  the  spectrum  of  a  complex  matrix  (see,  for  example, 
[5]).  In  this  paper,  the  QR  method,  as  described  in  [10],  is  used  for  computing  the 
eigenvalues  of  a  complex  linear  operator  on  a  hnite-dimensional  space. 


2.3.  Least  Squares  Newtou’s  Method.  Newton’s  method  is  a  well-known  numerical 
technique  for  solving  equations  of  the  form  F{x)  =  0,  where  F  :  is  a 

continuously  differentiable  function  of  the  form 


(2.1) 


F{x) 


^  /i(a^)^ 
\fnix)/ 


and  X  =  {xi, ,  XmY' ■  The  method  uses  the  Jacobian  matrix  J  of  F,  which  is  defined 
by  the  formula 


(2.2) 


.J{x) 


tw  ■ 

. .  TIl 

dXm 

dfn 

dXm 

{x)\ 

{x)J 


and  its  pseudoinverse,  which  is  denoted  by  J^{x).  We  illustrate  the  application  of  the 

method  in  the  following  theorem,  under  the  constraints  relevant  to  our  paper.  When 

m  ^  n,  the  procedure  is  referred  to  as  the  least  sguares  Newton’s  method  [11]. 
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Theorem  2.1.  Suppose  that  F  :  — >  M”  (m  >  n)  is  a  continuously  differentiable 

function  and  ^  G  satisfies  the  condition 

(2.3)  i"(O=0. 

Furthermore,  suppose  that  J{x)  defined  by  (2.2)  has  full  row  rank  for  all  x  in  an  e- 
neighborhood  ne{i)  =  {x  :  ||a;  —  .^||  <  e}  of  Let  {xk)k>Q  be  a  sequence  of  vectors 
defined  by  the  formula 

(2.4)  Xk+i  =  Xk-  j\xk)F{xk) , 
and  let  {hn)n>o  be  the  sequence  of  differences 

(2.5)  hn  Xji  ■ 

Then  there  exist  e  >  0,  5  >  0  and  a  nonneqative  integer  such  that,  if  ||a;o  —  ^||  <  e, 
then  the  inequality 

(2.6)  \\hn+l\\  <  S\\hn\\" 

is  true  for  all  k  >  N^^s- 

2.4.  Quadratures.  Quadrature  rules  are  expressions  of  the  form 

n 

(2.7)  '^Wiip{xi), 

i=l 

where  the  points  Xi  G  and  the  coefficients  tCj  G  M  are  referred  to  as  the  nodes  and 
weights  of  the  quadrature,  respectively.  They  serve  as  approximations  to  integrals  of 
the  form 

(2.8)  /  (p{x)w{x)  dx  , 

Jn 

where  w  is  a  positive  weight  function  and  C  is  the  region  of  integration. 

Quadratures  are  typically  chosen  so  that  the  formula  (2.7)  is  equal  to  the  desired 
integral  (2.8)  for  some  set  of  functions,  commonly  polynomials  of  some  hxed  order.  A 
quadrature  which  integrates  exactly  all  polynomials  of  (total)  degree  at  most  k  and 
does  not  integrate  exactly  some  polynomial  of  degree  /c  +  1  is  said  to  have  order  k.  For 
instance,  Gaussian  quadrature  rules  consist  of  n  nodes  and  have  order  2n  —  1. 

Given  a  weight  function  w,  an  integration  region  G  and  m  functions  dehned  on  G, 
a  classical  problem  in  numerical  integration  is  to  determine  the  minimum  number  of 
nodes  n  needed  to  integrate  exactly  all  the  m  functions,  and  to  hnd  quadratures  which 
achieve  that  minimum.  Such  quadratures  are  called  optimal.  In  one  dimension,  Gauss¬ 
ian  quadratures  are  known  to  be  optimal  rules  for  integrating  polynomials  up  to  a  hxed 
degree.  In  higher  dimensions,  however,  the  problem  has  not  been  solved  [2]. 

Since  each  node  a;  G  is  determined  by  d  parameters  and  one  additional  parameter 
is  given  by  the  weight  corresponding  to  the  node,  an  n-point  quadrature  in  is  deter¬ 
mined  by  n{d  +  1)  parameters.  For  this  reason,  it  is  natural  to  expect  that  an  optimal 
quadrature  in  d  dimensions  which  integrates  m  functions  has  n  =  \-^i\  nodes,  where 
[x]  denotes  the  smallest  integer  greater  or  equal  to  x.  Based  on  this  observation,  the 

authors  of  [11]  introduce  an  efficiency  coefficient  e  to  differentiate  quadratures  which 
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integrate  the  same  set  of  functions  and  have  a  different  number  of  nodes.  If  n  is  the 
number  of  nodes  of  the  quadrature,  m  is  the  number  of  functions  which  are  integrated 
exactly,  and  d  is  the  dimension  of  the  space,  then  the  coefficient  e  is  dehned  as 

(2.9)  e  =  "" 


n{d  +  1) 


For  Gaussian  quadratures  in  one  dimension,  e  is  equal  to  1.  Quadratures  having  the 
coefficient  e  close  to  1  are  called  efficient. 


2.5.  Tensor  Product  Rules.  An  obvious  way  of  obtaining  quadratures  in  two  dimen¬ 
sions  is  to  combine  two  one-dimensional  quadratures  in  the  following  way. 

Let  R  be  an  m-point  quadrature  over  an  interval  /, 

m  p 

(2.10)  R{f)  =  ^Wif{xi)^  f{x)dx, 

i=i  ^ 

and  S  be  an  n-point  quadrature  over  an  interval  J, 

m  p 

(2.11)  S{f)  =  ^Vjf{xj)^  f{x)dx. 

j=i  JJ 

The  tensor  product  of  R  and  S  is  the  mn-point  quadrature  rule  over  the  rectangle  I  x  J 
dehned  by 

m  n 

(2.12)  R  X  S{f)  =  EE 

i=l  j=l 

The  following  theorem  follows  from  Fubini’s  Theorem;  a  proof  can  be  found,  for  instance, 
in  [4]. 


f{x,y)  dxdy. 


'Ixj 


Theoreru  2.2.  Let  R  be  a  quadrature  rule  that  integrates  f{x)  exactly  over  the  interval  I 
and  S  be  a  quadrature  rule  that  integrates  g{y)  exactly  over  the  interval  J .  Furthermore, 
let  h{x,y)  =  f{x)g{y).  Then  the  tensor  product  quadrature  R  x  S  integrates  h{x,y) 
exactly  over  I  x  J. 


Tensor  product  quadratures  are  relatively  easy  to  construct,  however,  they  are  not 
optimal  in  terms  of  efficiency.  For  example,  let  us  consider  two  n-point  Gaussian  quadra¬ 
tures  on  the  interval  [—1, 1],  each  of  them  of  order  2n  —  1.  Their  tensor  product  is  an 
n^-point  quadrature  on  the  square  [—1, 1]  x  [—1,1]  of  order  2n  —  1  (all  monomials  x^y^ 
of  total  degree  i  +  j  at  most  2n  —  1  can  be  expressed  as  products  of  monomials  in  the 
variables  x  and  y  of  degree  at  most  2n  —  1,  hence  they  are  integrated  exactly  by  the 
tensor  product;  however,  there  are  polynomials  of  degree  2n,  for  instance  x'^^,  which 
are  not  integrated  exactly  by  this  quadrature).  The  number  of  monomials  x^y^  of  total 
degree  at  most  2n  —  1  is  n(2n  -|-  1),  therefore  the  efficiency  coefficient  of  the  tensor 
product  of  two  Gaussian  quadratures  is 

n{2n  -1-1)  2n  -I-  1 
3n 

quantity  which  approaches  |  as  n  increases. 
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3.  Mathematical  Apparatus 


The  main  result  of  this  section  is  Theorem  3.2,  which  allows  us  to  discretize  bounded 
convex  regions  in  the  complex  plane.  It  is  based  on  the  observation  that  all  the  eigen¬ 
values  of  a  certain  linear  operator  on  a  hnite-dimensional  subspace  of  L‘^{X),  where  X 
is  a  convex,  closed  and  bounded  set  in  the  plane,  fall  within  X.  Our  proof  of  Theorem 
3.2  involves  the  geometrical  result  introduced  in  the  following  section. 

3.1.  A  geometrical  lemma. 

Lemma  3.1.  Let  X  be  a  convex,  closed  and  bounded  subset  of  the  complex  plane  and 
\  E  C  be  a  point  outside  X .  Then  there  exists  a  point  p  G  C  such  that  pi  lies  outside  X 
and 

(3.1)  1-^  ~  hi  >  k  ~  hi 
for  all  z  in  X . 

Proof.  Let  us  denote  by  d  the  distance  from  A  to  the  set  X: 

(3.2)  d  =  d{X,X)  =  mi\z-X\. 

Since  X  is  closed,  the  inhmum  d  is  attained  for  some  ;2a  in  X,  and  is  positive  because 
A  is  not  in  X: 

(3.3)  d=\zx-X\>0. 

Consider  the  circle  C  of  radius  d  centered  at  A  and  the  tangent  /  to  C  at  the  point 
2:a  (see  Figure  1). 


Figure  1.  The  tangent  I  to  the  circle  C  centered  at  X,  at  the  point  z\. 
The  point  zx  minimizes  the  distance  from  X  to  the  set  X . 
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The  line  /  determines  two  half-planes:  one  containing  the  point  A  and  one  not  contain¬ 
ing  the  point  A.  We  observe  that  X  has  no  points  inside  the  open  half-plane  determined 
by  /  that  contains  the  point  A. 

Let  us  prove  this  observation  by  contradiction.  To  this  end,  suppose  that  there  exists 
a  point  z  in  X  which  lies  in  the  open  half-plane  containing  A.  Since  2:  and  A  lie  in  the 
same  half-plane  determined  by  the  tangent  I,  and  2:  is  not  on  the  line  /,  a  part  of  the 
segment  [zz\]  falls  inside  the  circle  C  (see  Figure  2). 


Figure  2.  Proof  of  the  observation  by  contradiction:  if  a  point  z  in  X 
lies  in  the  open  half-plane  determined  by  the  line  I  which  contains  the 
point  \,  then  some  points  of  X  fall  inside  the  circle  C . 

However,  X  is  convex  and  the  points  2:  and  2:a  lie  in  X,  therefore  the  whole  line 
segment  [zzx]  is  contained  in  X.  Consequently,  some  points  of  X  fall  inside  the  circle 
C,  at  distance  less  than  d  to  A,  contradicting  our  choice  of  2:a.  We  conclude  that  the 
set  X  is  contained  in  the  closed  half-plane  determined  by  /  that  does  not  contain  A. 

Since  X  is  bounded,  we  can  hnd  a  square  box  L1L2T3L4  that  contains  X.  Moreover, 
by  the  previous  observation,  we  can  choose  the  box  such  that  the  points  Li,  L2  lie  on  the 
line  /  and  the  points  L3,  L4  he  in  the  half-plane  determined  by  I  that  does  not  contain 
A.  Finally,  we  can  choose  the  square  box  big  enough  so  that  the  length  s  of  one  of  its 
sides  is  greater  than  d  (see  Figure  3). 

Let  fi  be  the  point  on  the  the  half-line  passing  through  the  points  A  and  2a  with 
initial  point  A  with  the  following  property: 

(3.4)  \zx-p\  =  ^. 

2 

Since,  by  construction,  s  >  d,  multiplying  by  |  we  obtain  ^  >  s,  hence  the  point  p 

lies  outside  the  square  box  L1L2L3L4  and  therefore  outside  X  (see  Figure  4). 
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Figure  3.  X  is  enclosed  in  a  square  box  L1L2F3L4  lying  in  the  half¬ 
plane  determined  by  1  that  does  not  contain  A.  The  length  s  of  a  side  of 
the  square  is  greater  than  d. 

Li 

1 

L4 

.  d  \ 

^  ) 

k  / 

F 

d 

L2 

L3 

Figure  4.  We  define  pi,  to  be  the  point  on  the  half-line  passing  through 

2 

the  points  X  and  z\  with  initial  point  A  such  that  pi  is  at  distance  ^  from 
2 

2:a.  Since  ^  >  s,  pi  lies  outside  the  box  L1L2F3L4  and  therefore  outside 

Let  us  show  that  /i  also  satishes  the  other  condition  in  the  hypothesis,  i.e., 


(3.5)  1-^  ~  hi  >  k  ~  hi  ) 

for  all  z  in  X. 

Indeed,  let  2:  be  an  arbitrary  point  in  X  and  let  be  its  projection  on  the  line  passing 
through  the  points  and  p.  The  distance  from  z  to  fi  is  the  length  of  the  hypotenuse 
in  the  right  triangle  Azz^n  (see  Figure  5). 


Figure  5.  The  distance  from  a  point  z  in  X  to  p,  is  the  length  of  the 
hypotenuse  in  the  right  triangle  Azz^p,  where  z^  is  the  projection  of  z  on 
the  line  passing  through  the  points  zx  and  p 


Since  p  lies  outside  the  square  L1L2L3L4,  the  point  2;^  lies  on  the  line  segment  [zxp], 
hence 

(3.6)  \z^  -  hi  <  I^A  -  hi  • 

Note  that  since  X  is  contained  in  the  square  L1L2L3L4  and  the  line  zxp  passing 
through  the  points  2:a  and  p  is  perpendicular  on  the  side  L1L2,  the  projection  Zf^  of  the 
point  2:  in  X  on  zxp  falls  inside  L1L2L3L4.  Since  the  largest  distance  between  any  two 
points  inside  a  square  is  the  length  of  the  diagonal,  we  have 

(3.7)  \Zfj^  —  z\  <  sa/2  . 

By  the  Pythagorean  Theorem  and  inequalities  (3.6)  and  (3.7), 

(3.8)  |;2  -  /i|^  =  \Zf,  -  /ip  +  \Zf,  -  zj^  <  jzx  -  hP  +  . 
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On  the  other  hand,  by  the  equation  dehning  /i  (3.4)  we  have  the  following  formula 
for  the  square  of  the  length  of  the  line  segment  [A/i]: 

(3.9)  |A  —  /i|^  =  (|A  —  z\\  +  I^A  ~  =  {d  +  \zx  —  /i|)^ 

(3.10)  =  +  2d\zx  —  Ail  +  \zx  —  /iP  =  +  2d—  +  \zx  —  /iP 

d 

(3.11)  =  d‘^  +  2s^  +  \zx  -  . 

Since  d  is  positive,  by  formula  (3.11)  we  have 

(3.12)  |A  -  A^r  =  d^  +  2s^  +  |;2a  -  >  2s^  +  |;2a  -  A^r  • 

Combining  inequalities  (3.12)  and  (3.8)  we  obtain 

(3.13)  |A-/ip  > 

Finally,  taking  the  square  root  of  both  sides  yields 

(3.14)  |A  —  Ai|  >  k  —  Ai|  • 

Since  z  in  X  was  chosen  arbitrarily,  inequality  (3.14)  holds  for  all  in  X,  therefore 
/i  is  a  point  satisfying  the  conditions  in  the  hypothesis. 

□ 

3.2.  Analytical  Theorem.  In  the  beginning  of  this  section,  we  summarize  a  few  well- 
known  facts  in  functional  analysis,  which  are  used  in  the  proof  of  Theorem  3.2. 

Let  X  be  a  measurable  subset  of  the  complex  plane.  The  linear  space  of  square 
integrable  complex  functions  on  X  is  generally  denoted  by  L^{X).  As  is  well-known, 
L‘^{X)  has  the  structure  of  a  complex  Hilbert  space  with  the  inner  product  dehned  by 


(3.15) 


(f,g)  ^  f(x  ,y)g{x,y)  dxdy^ 


Let  U  be  an  arbitrary  finite-dimensional  subspace  of  L‘^{X).  We  denote  by  Pu  the 
orthogonal  projection  L^{X)  — U  and  by  Tu  the  operator  multiplication  hj  x  +  iy  on 
the  space  U : 


(3.16a)  Tu:  U — ^  ^^(X) , 

(3.16b)  Tu{f{x,  y))  =  {x  +  iy)  ■  f{x,  y) . 


It  is  easy  to  see  that  Tu  is  a  linear  operator  on  A/,  hence  the  composition  Pu  o  T[/  is  a 
linear  operator  on  the  hnite-dimensional  space  U. 

The  following  theorem  lists  a  few  basic  properties  of  the  orthogonal  projection  oper¬ 
ator.  A  proof  can  be  be  found,  for  instance,  in  [1]. 


Theorem  3.1.  Let  X  be  a  measurable  subset  ofC  and  U  be  a  nonzero  finite- dimensional 
subspaee  of  Lfi{X).  The  orthogonal  projection  Pu  has  the  following  properties: 

{1)  Pu  is  a  linear  operator  with  operator  norm  equal  to  1. 

(2)  The  restriction  of  Pu  to  the  space  U  is  the  identity  operator  on  U ,  i.e.,  Pu{t)  = 
(f  for  all  functions  if  in  U . 
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(3) 

(3.17) 


The  operator  Pu  is  self-adjoint,  i.e., 


for  all  in  Lf{X), 
{p>,Pu{f)))  . 


The  principal  tool  for  the  constructions  of  discretizations  and  quadratures  in  this 
paper  is  the  following  analytical  theorem. 

Theorem  3.2.  Let  X  be  a  measurable,  convex,  closed  and  bounded  subset  ofC,  and  U 
be  a  finite- dimensional  linear  subspace  of  L‘^{X).  Let  Tjj  be  the  linear  operator  defined 
by  (3.16),  and  Pu  be  the  projection  onto  the  space  U.  Then  all  the  eigenvalues  of  the 
operator  Pu  o  Tu  fall  inside  X . 

Proof.  The  theorem  is  proved  by  contradiction.  To  this  end,  suppose  that  the  operator 
Pu  o  Tu  has  an  eigenvalue  A  G  C  outside  X.  By  Lemma  3.1,  there  exists  /i  G  C  such 
that  l/i  —  A|  >  \fi  —  z\,  for  all  z  E  X,  hence 

(3.18)  Ih  ~  -^1  >  iiiax  In  —  z\ . 

zSX 

Since  A  is  an  eigenvalue  for  Pu  oTu,  there  exists  a  nonzero  (p  E  U  such  that 

(3.19)  Pu{Tu{(p))  =  X(p 
Subtracting  pep  from  both  sides  of  (3.19)  we  obtain 

(3.20)  Pu{Tu{(p))  -  pap  =  X(p  -  pop  . 

Since  (p  belongs  to  the  subspace  U,  by  Theorem  3.1,  we  have  p,ip  =  Pu{pp>),  therefore 
we  can  rewrite  (3.20)  as 

(3.21)  PuiTuip))  -  PuiTP)  =  X(p  -  p,(p  , 
which  is  equivalent  to 

(3.22)  Pu  (Tuip)  -  p.p)  =  (X  -  p,)p  . 

The  operator  norm  of  the  projection  Pu  is  1  by  Theorem  3.1,  therefore 

(3.23)  \\Tu{p)  -  TtW  >  \\Pu{Tu{p)  -  pp)\\  . 

Combining  equation  (3.22)  with  inequality  (3.23)  we  obtain 

(3.24)  \\Tuip)  -  l^p\\  >  II  (A  -  p)p\\  . 

Using  the  formula  dehning  Tu  (3.16b),  inequality  (3.24)  becomes 

(3.25)  \\{x  +  iy)p{x,y)  -  pp{x,y)\\  >  \\{X  -  p,)p\\  , 
which  is  equivalent  to 

(3.26)  \\{x  +  iy-  p,)p{x,y)\\  >  |A-/i|  •  \\pix,y)\\  . 

The  complex  variable  x  +  iy  takes  values  only  within  the  set  X,  therefore 

^-hlV  \\T{.x,y)\\  >  \\{x  +  iy  -  p)  p{x,y)\\ . 


(3.27) 


( 


max 
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Inequalities  (3.26)  and  (3.27)  together  yield 

(3.28)  ■  \\^{x,y)\\  >\\-y\-  \\(p{x,y)\\. 

However,  (p  is  nonzero,  hence  we  can  divide  both  sides  of  (3.28)  by  the  positive  number 
1199(0;,  y)  II  and  obtain 

(3.29)  max|;2  —  /i|  >  |A  —  /i| , 

z£X 

inequality  which  contradicts  our  choice  of  y  (3.18). 

□ 

3.3.  Remarks. 

Remark  3.1.  The  eigenvalues  of  the  operator  Pu  o  Tjj  are  complex  numbers  that 
fall  inside  X,  hence  they  give  us  a  discretization  of  the  region  X  that  encapsulates 
information  about  the  hnite-dimensional  space  U.  Theorem  3.2  provides  a  method  for 
generating  discretizations  of  (arbitrary)  bounded  convex  regions  in  two  dimensions.  For 
a  given  region,  one  can  vary  the  choice  of  the  subspace  U  in  order  to  obtain  different 
discretizations. 

Remark  3.2.  Let  G  be  the  group  of  symmetries  of  X,  i.e.,  G  consists  of  all  isometries 
under  which  X  is  invariant.  Suppose  that  the  subspace  U  of  L‘^{X)  is  invariant  under 
the  action  of  the  group  G,  i.e., 

(3.30)  g{f)^U, 

for  all  g  &  G  and  f  E  U,  where  the  action  of  an  isometry  in  G  on  a  function  f  in  U 
is  dehned  by 

(3.31)  g{f)  =  fgE  L\X) , 

(3.32)  fg{x,y)  =  f{g~\x,y)) ,  for  all  (x,  I/)  e  X  . 

Then  the  spectrum  of  the  operator  P  o  T  is  also  invariant  under  the  action  of  G,  i.e.,  if 
A  is  an  eigenvalue  of  P  o  T,  then  g{X)  is  also  an  eigenvalue  of  P  o  T. 

We  can  summarize  this  remark  in  the  following  way.  If  the  subspace  U  is  chosen  as  to 
be  invariant  under  the  action  of  the  symmetry  group  of  X,  then  the  discretization  given 
by  Theorem  3.2  will  be  symmetric  (invariant  under  the  symmetry  group  of  X). 

Remark  3.3.  The  nodes  of  the  discretization  given  by  Theorem  3.2,  together  with 
suitably  chosen  weights,  provide  a  good  initial  point  for  Newton’s  method,  which  can 
be  used  to  obtain  a  quadrature  rule  on  X.  This  application  of  Theorem  3.2  is  illustrated 
in  Section  4. 

Remark  3.4.  Let  {(pi,ip2,  ■  ■  ■  ,^n}  be  an  orthonormal  basis  of  U.  The  matrix  A 
representing  Pu  o  Tu  in  this  basis  has  entries  (0^^)  dehned  by 

(3.33)  ttij  =  ( (Pf/  o  Tu){(Pj),  Pi )  . 

By  Theorem  3.1,  the  orthogonal  projection  operator  is  self-adjoint,  hence 

(3.34)  ttij  =  {Tu{pj),  PuiPi))  ■ 
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Since  ^pi  belongs  to  U,  by  Theorem  3.1,  Pjj{ipi)  =  ipi,  for  all  1  <  i  <  n.  Therefore  we 
can  rewrite  equation  (3.34)  as 


(3.35) 


Finally,  we  can  use  the  dehnition  of  Tu  (3.16b)  and  of  the  inner  product  (3.15)  to  write 
equation  (3.35)  in  the  following  explicit  form; 


(3.36) 

(3.37) 

(3.38) 


=  {{x  +  iy)p>j{x,y),  (pi{x,y)) 


{x  +  iy)  -(fjix^y)  ■  (pi{x,y)  dxdy. 


X 


This  explicit  formula  is  used  in  Stage  3  of  the  algorithm  described  in  the  next  section. 


4.  Numerical  Algorithm 

This  section  describes  a  numerical  algorithm  for  the  construction  of  a  class  of  inter¬ 
polation  and  quadrature  formulae  on  regular  polygons  in  the  plane.  The  algorithm  can 
be  generalized  to  arbitrary  bounded  convex  regions  in 

For  simplicity,  we  assume  that  the  regular  polygon  is  centered  at  the  origin  and  its 
vertices  lie  on  the  unit  circle.  Let  s  be  the  number  of  sides  of  the  polygon  (s  >  3),  and 
X  be  the  convex  region  enclosed  by  the  polygon.  For  a  given  positive  integer  n,  we 
denote  by  the  space  of  polynomials  on  X  of  total  degree  at  most  n: 

(4.1)  “Pn  =  {P{x,y)  :  total  degree  of  P  <  n}  C  L^{X) . 

One  can  easily  see  that  the  dimension  of  Pn  is  given  by  the  formula 

(4.2)  dim!P„  =  (n -I- l)(n -I- 2)/2  . 

The  algorithm  uses  Theorem  3.2  with  U  =  Pn,  P  '■  LP{X)  — >  U  the  orthogonal 
projection  operator,  and  T  the  multiplication  by  [x+iy)  operator  (as  dehned  by  (3.16b)), 
in  order  to  discretize  the  region  X.  Note  that  the  space  Pn  is  invariant  under  the  action 
of  the  dihedral  group  Dg,  the  symmetry  group  of  X.  Therefore,  by  Remark  3.2,  the 
discretization  given  by  the  theorem  is  symmetric  (invariant  under  the  action  of  Dg). 
The  nodes  of  this  discretization  are  used  as  initial  nodes  for  Newton’s  method  in  order 
to  construct  quadrature  rules  on  the  regular  polygon  X. 

The  algorithm’s  input  is 

s  =  the  number  of  sides  of  the  polygon 
n  =  the  initial  maximum  degree  of  the  polynomials 
m  =  the  maximum  degree  of  the  polynomials  to  be  integrated. 

By  equation  (4.2),  the  dimensions  of  Pn  and  Pm  are  N  =  (n  -|-  l)(n  -|-  2)/2  and  M  = 

(m  -|-  l)(m  -|-  2)/2,  respectively.  The  algorithm’s  output  is  a  quadrature  rule  on  X 
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consisting  of  a  set  of  N  nodes  (a;i,  |/i),  (0:2, 1/2),  •  •  • ,  (a^TV,  l/iv)  and  a  set  of  N  weights 
Wi,W2, . . . ,  wj\f  snch  that 


N 

(p{x,  y)  dx  dy  ~  E  (p{xj,yj)wj 

t=i 


for  any  polynomial  ip  in  the  M-dimensional  space  Tm- 


The  algorithm  proceeds  in  4  stages. 

In  the  hrst  stage,  we  hnd  a  quadrature  of  the  polygon  by  dividing  it  into  triangles,  and 
then  constructing  a  tensor  product  of  Gaussian  quadratures  on  each  of  the  triangles.  In 
the  second  stage  we  construct  an  orthonormal  basis  for  the  space  iteratively,  using 
the  SVD.  In  the  third  stage  we  hnd  the  eigenvalues  of  the  operator  PoT.  In  the  fourth 
stage  we  use  Newton’s  method  to  modify  the  position  of  the  eigenvalues  and  of  the 
initial  weights,  found  using  the  Least  Squares  Method,  in  order  to  obtain  a  quadrature 
rule  which  integrates  all  polynomials  of  total  degree  at  most  m. 


Stage  1:  Discretize  the  polygon. 

In  this  stage  we  hnd  a  quadrature  of  the  regular  polygon  by  performing  the  following 
sequence  of  steps: 

(1)  Divide  the  regular  polygon  A1A2  . . .  An  centered  at  O  into  n  triangles 
OA1A2,  OA2A3, . . . ,  OAnAi,  where  O  is  the  origin  of  the  axes. 

(2)  Discretize  each  of  the  triangles  by  constructing  the  tensor  product  of  two  one¬ 
dimensional  Gaussian  quadratures  on  an  interval.  The  order  of  the  Gaussian 
quadratures  is  chosen  sufficiently  large  so  that  the  tensor  product  quadrature 
integrates  all  polynomials  of  total  degree  at  most  m  exactly. 

(3)  Assemble  all  the  nodes  and  weights  from  all  the  triangles  to  obtain  a  quadrature 
of  the  polygon. 

Stage  2:  Construct  the  orthonormal  basis. 

In  this  stage,  we  hnd  an  orthonormal  basis  for  the  space  The  basis  is  used  for 
determining  the  eigenvalues  of  the  operator  =  PoT  on  the  space  therefore 

the  choice  of  this  basis  is  not  important  -  any  basis  yields  the  same  eigenvalues.  The 
following  algorithm  for  constructing  an  orthonormal  basis  is  not  optimal,  however,  it  is 
stable  and  easy  to  implement. 

(1)  Gonsider  the  vector  corresponding  to  the  constant  polynomial  1.  Normalize  it. 
The  resulting  vector  provides  an  orthonormal  basis  for  Tq- 

(2)  Suppose  we  have  constructed  an  orthonormal  basis  of  Tk,  consisting  of  K  = 
{k  +  1)(A;  -|-  2)/2  vectors.  Multiply  each  of  the  vectors  in  this  basis  by  x  and 
y,  respectively,  and  put  the  resulting  vectors,  together  with  the  vectors  in  the 
basis  of  Tfc,  in  a  matrix  having  3K  columns.  By  performing  SVD  on  this  matrix 
we  hnd  an  orthonormal  basis  for  its  column  space,  which  is  also  an  orthonormal 
basis  for  Tfc+i- 

(3)  Apply  iteratively  the  procedure  described  in  the  previous  step  to  obtain  an 
orthonormal  basis  for 
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Stage  3:  Find  the  eigenvalues  of  P  o  T. 

(1)  Compute  the  entries  {aij)  of  the  matrix  A  representing  P  o  T  in  the  basis 
{'•Pi)  P2)  ■  ■  ■ )  Pn}  constructed  in  Stage  2  using  Formula  (3.38): 


aij  — 


{x  +  iy)pi{x,  y)pj{x,  y)  dx  dy . 


X 

The  integrals  are  evaluated  using  the  quadrature  constructed  in  Stage  1. 

(2)  Find  the  eigenvalues  Ai,  A2, . . . ,  Atv  of  the  complex  matrix  A  using  the  QR  algo¬ 
rithm  described  in  Section  2.2. 


Stage  4:  Construct  the  quadrature. 

In  this  stage,  we  hnd  a  quadrature  rule  consisting  of  N  nodes  and  N  weights  which 
integrates  all  M  polynomials  of  total  degree  at  most  m.  The  algorithm  uses  Newton’s 
method  and  consists  of  the  following  seqnence  of  steps: 

(1)  Constrnct  an  orthonormal  basis  (pi,(p2, . . .  ,Pm  of  the  space  using  the  pro- 
cednre  described  in  Stage  2. 

(2)  The  initial  qnadrature  points  {x^,  y^),  1/2),  •  •  • ,  {x%)  y^)  for  Newton’s  method 

are  given  by  the  eigenvalues: 


+  =  J  =  1,2,  ...,iV. 

The  initial  quadrature  weights  10^,102,  ■■■ ,  w%  are  found  by  solving  the  following 
linear  system  using  Least  Squares: 


Pji^yi)^!  +  Pj{xl,y2)w2  +  •  •  •  +  Pj{x%,y%)w%  =  Ij, 


J  =  l,2, 


M, 


where 

=  Jj  Pjix,y)dxdy,  j  =  l,2,...,M 

are  the  exact  values  of  the  integrals  of  the  polynomials  in  the  orthonormal  basis 
of  Tm  fonnd  in  step  (1).  The  integrals  Ij,j  =  1,2, ...  ,M  are  fonnd  using  the 
initial  quadratnre  constrncted  in  Stage  1. 

(3)  Find  a  quadratnre  rule  for  integrating  all  polynomials  in  This  is  equivalent 
to  solving  the  following  nonlinear  system  with  3N  nnknowns 


{wi,W2,  .  .  .,Wn,Xi,X2,  .  .  .  ,XN,yi,y2,  •  •  •  ,  |/v) 
and  M  equations: 


Pj{xi,yi)wi  +  Pj{x2,y2)w2  -f  .  .  .  +  Pj{xN,yN)wN  =  Ij)  j  =  1,2, . . .  ,M. 


Finding  a  solntion  to  this  system  is  eqnivalent  to  hnding  a  zero  of  the  smooth 
fnnction  F  :  defined  by 


F{wi, . . . ,  WN)  a^i, . . . ,  XN)  yi)---)  yx) 


/  N  \ 

I  y^^pi{xi,yi)wi  -  h 

i=i 


N 


'y  ^  Piiiixi)  yi 

\i=l 


\Wi 


l-M 
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In  order  to  find  a  zero  of  F,  we  use  Least  Squares  Newton’s  method  (as  described 
in  Section  2.3)  with  initial  point 

(w?,  w%,  x?,  xl,...,  x%,  vl,...,  y%) . 

The  weights  and  nodes  given  by  this  zero  of  F  dehne  an  iV-point  quadrature 
which  integrates  all  polynomials  in  the  M-dimensional  space  Tm- 

5.  Numerical  Results 

We  implemented  the  procedure  described  in  Section  4  and  tested  it  on  a  number  of  ex¬ 
amples.  Several  of  the  results  are  listed  below,  in  order  to  demonstrate  the  performance 
of  the  algorithm.  The  algorithm  was  implemented  in  FORTRAN  and  compiled  with  the 
Lahey-Fujitsu  FORTRAN  95  compiler.  Computations  were  performed  in  double  and 
extended  precision  arithmetic. 

5.1.  Quadratures  and  interpolation  formulae  on  the  triangle.  In  this  section  we 
describe  the  quadratures  and  interpolation  formulae  on  an  equilateral  triangle  generated 
by  the  algorithm.  The  triangle  is  centered  at  the  origin  and  has  vertices  located  at  the 
points  {1,  on  the  unit  circle.  The  following  procedure  is  used. 

We  begin  with  an  initial  polynomial  degree  n.  Using  the  algorithm  described  in 
Section  4,  we  find  the  maximum  polynomial  degree  m  for  which  Least  Squares  Newton’s 
method  yields  a  quadrature  having  positive  weights.  The  resulting  quadratures  have 
N  =  (n  -|-  l)(n  -|-  2)/2  nodes  and  can  integrate  all  the  polynomials  in  the  space  of 
dimension  M  =  [m  +  l){m  +  ‘2)/2. 

Remark  5.1.  Given  the  orthonormal  basis  991,  (/)2,  •  •  • ,  of  the  space  Tn  constructed 
in  Stage  2  of  the  algorithm,  the  N  nodes  Xi,X2,  ■  ■  ■  ,xn  of  the  quadrature  yield  an 
interpolation  formula  relative  to  the  functions  {(pi\.  The  number  of  nodes  for  this 
interpolation  formula  equals  the  number  of  functions.  The  formula  is  stable  if  the 
condition  number  of  the  associated  interpolation  matrix  B  is  small,  where  B  is  dehned 
by  the  equations 

(5.1)  B  =  {hij),  bij  =  (pi{xj),  i,  j  =  1,2, . . . ,  W 

In  Tables  1  and  2  we  list  the  quadratures  generated  by  this  procedure,  together  with 
the  corresponding  efficiency  coefficients  e  (as  dehned  by  Equation  2.9)  and  condition 
numbers  of  the  interpolation  matrices  B  (as  dehned  by  Equation  5.1).  The  quadratures 
of  degrees  9,  11,  16,  21,  26,  29  and  31  have  rotational  symmetry,  but  not  rehection 
symmetry,  while  all  the  remaining  quadratures  are  fully  symmetric  —  cf.  Remark  5.2. 

Remark  5.2.  In  most  cases,  since  the  eigenvalues  are  invariant  under  the  action  of  the 

symmetry  group  of  the  triangle,  modifying  their  position  by  Least  Squares  Newton’s 

method  yields  quadratures  on  the  triangle  which  are  fully  symmetric.  In  a  few  instances, 

however,  relaxing  the  requirement  that  the  quadratures  have  full  symmetry,  and  allowing 

only  rotational  symmetry  instead,  allows  us  to  obtain  quadratures  for  higher  orders  m. 

The  phenomenon  has  not  been  investigated  further. 
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n 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

m 

2 

4 

5 

7 

9 

11 

12 

14 

16 

17 

N 

3 

6 

10 

15 

21 

28 

36 

45 

55 

66 

M 

6 

15 

21 

36 

55 

78 

91 

120 

153 

171 

e 

0.67 

0.83 

0.70 

0.80 

0.87 

0.93 

0.84 

0.89 

0.93 

0.86 

CN 

1.0 

1.0 

1.2 

1.4 

1.6 

2.1 

1.8 

2.3 

4.2 

5.9 

Table  1.  Triangle  quadratures  -  part  1.  For  each  initial  polynomial 
degree  n  from  1  to  10,  we  list  the  corresponding  maximum  degree  m 
of  a  quadrature  having  positive  weights.  The  number  of  nodes  of  the 
quadrature  is  iV  =  (n  +  l)(n  +  2)/2  and  the  number  of  polynomials  which 
can  be  integrated  exactly  is  M  =  dimT^  =  (m  +  l)(m  +  2)/2.  The 
efficiency  of  the  quadrature  is  e  =  M/{3N).  The  condition  number  of  the 
interpolation  matrix  (5.1)  is  listed  in  the  row  labeled  CN. 


n 

11 

12 

13 

14 

15 

16 

17 

18 

19 

m 

19 

21 

22 

23 

26 

27 

29 

31 

32 

N 

78 

91 

105 

120 

136 

153 

171 

190 

210 

M 

210 

253 

276 

300 

378 

406 

465 

528 

561 

e 

0.90 

0.93 

0.88 

0.83 

0.93 

0.88 

0.91 

0.93 

0.89 

CN 

5.0 

9.7 

6.3 

13.1 

68.6 

144.0 

67.8 

71.5 

56.2 

Table  2.  Triangle  quadratures  ~  part  2.  For  each  initial  polynomial 
degree  n  from  11  to  19,  we  list  the  corresponding  maximum  degree  m 
of  a  quadrature  having  positive  weights.  The  number  of  nodes  of  the 
quadrature  is  =  (n  +  l)(n  +  2)/2  and  the  number  of  polynomials  which 
can  be  integrated  exactly  is  M  =  dimT^  =  {m  +  l){m  +  ‘2)12.  The 
efficiency  of  the  quadrature  is  e  =  M/{3N).  The  condition  number  of  the 
interpolation  matrix  (5.1)  is  listed  in  the  row  labeled  CN. 


In  Tables  3  and  4,  we  list  two  of  the  quadrature  rules  on  the  triangle  obtained  using 
the  procedure  described  above.  Since  each  of  the  quadratures  has  full  symmetry,  only  a 
subset  of  generating  nodes  and  corresponding  weights  is  sufficient  for  describing  each  of 
the  quadratures  (see,  for  example,  [11]  for  more  details).  The  remaining  nodes  may  be 
obtained  via  the  action  of  the  symmetry  group  of  the  triangle,  with  the  obvious  obser¬ 
vation  that  points  belonging  to  the  same  orbit  have  the  same  weight.  The  quadratures 
were  computed  in  extended  precision  and  are  listed  with  20  digits. 

In  Figure  6,  we  show  a  picture  of  the  nodes  of  the  fully  symmetric  quadrature  of  order 
32. 

5.2.  Quadratures  on  other  polygons.  Constructing  quadratures  on  regular  poly¬ 
gons  having  at  least  4  sides  using  the  algorithm  described  in  Section  4  is  slightly  more 
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Weight 


-.16868265317233866793E+00 

-.15808488336403957216E+00 

-.35347157488571107334E+00 

-.47009724117048709338E+00 

-.36136906345529792533E+00 

-.47259519998535818893E+00 

-.47234427519430125315E+00 

-.35112713454811437607E+00 

-.47085943420511108264E+00 

-.47949927002953977850E+00 


0.29216692565001003106E+00 

-.12519480877035419803E-30 

0.14203631281845219952E+00 

0.41794071711619191254E-30 

0.40253738576591622170E+00 

0.24682585755277199100E+00 

0.50374523862970667221E+00 

0.60817003695340733938E+00 

0.70900463999219136033E+00 

0.83051709788335153417E+00 


0.58935620335540765011E-01 

0.67452219369552326397E-01 

0.44503540566573551812E-01 

0.17897106324896028214E-01 

0.36144179089957608755E-01 

0.18293541893943411802E-01 

0.16917794610623270344E-01 

0.28921171272990801675E-01 

0.12279549423359387150E-01 

0.35293734203249423613E-02 


Table  3.  A  fully  symmetric  45-point  quadrature  of  order  14  on  the  tri¬ 
angle.  Only  10  generating  nodes  are  listed;  the  remaining  nodes  may  be 
obtained  via  the  action  of  the  group  of  symmetries  of  the  triangle. 


Weight 


-.12399257021092653653E+00 

-.13177700392715291090E+00 

-.28778371647987152368E+00 

-.40630672781840910415E+00 

-.29724898552695112330E+00 

-.41356124787667862272E+00 

-.48264071535159024307E+00 

-.41408886151213149256E+00 

-.27842964300896413291E+00 

-.48370152343347647465E+00 

-.48281402069237875749E+00 

-.40528714204595318336E+00 

-.48105533677985120342E+00 

-.41966365447253047169E+00 

-.48630312842497788637E+00 

-.48251983199382645129E+00 


0. 19374825509442976738E-28 
0.22824446607103231387E+00 
0.10998294723245009451E+00 
-.14456488256121493556E-28 
0.31498865266126179066E+00 
0.16719641141490991682E+00 
0.95221898752783636336E-01 
0.38336091989130268782E+00 
0.48225428802479053639E+00 
0.30142300793320924164E+00 
0.49942229259718307221E+00 
0.57394856816236337313E+00 
0.64968083098649002499E+00 
0.72687877163645268678E+00 
0.76255697108968090949E+00 
0.83574886467290610008E+00 


0.41450918407684045121E-01 
0.38440468676165286217E-01 
0.31317527185910258616E-01 
0.14002753888728943368E-01 
0.26827201883108899153E-01 
0.19961899945473734702E-01 
0.86539518883907294556E-02 
0.20240988214359192248E-01 
0.23989433658357479662E-01 
0.87991644027906763838E-02 
0.78178774263925103734E-02 
0.18008780871442209707E-01 
0.60973338938233607918E-02 
0. 105471 16649857438346E-01 
0.34615238341665281813E-02 
0.22095115197099314446E-02 


Table  4.  A  fully  symmetric  78-point  quadrature  of  order  19  on  the  tri¬ 
angle.  Only  16  generating  nodes  are  listed;  the  remaining  nodes  may  be 
obtained  via  the  action  of  the  group  of  symmetries  of  the  triangle. 


involved.  In  contrast  to  the  triangle  case,  the  operator  P  o  T  has  sometimes  multiple 
eigenvalues;  in  all  cases,  the  only  multiple  eigenvalue  is  0,  the  center  of  the  polygon. 

It  is  easily  seen  that,  in  some  cases,  0  is  a  multiple  eigenvalue  of  the  operator  P  o  T 
via  the  following  argument. 

Observation  5.1.  Suppose  that  the  number  N  of  eigenvalues  is  not  divisible  by  the 
number  of  sides  s.  By  Remark  3.2,  the  set  of  eigenvalues  is  invariant  under  rotations, 
hence  we  can  divide  the  nonzero  eigenvalues  into  s  groups,  each  corresponding  to  one 
of  the  sides  of  the  polygon.  Consequently,  there  are  s  ■  k  nonzero  eigenvalnes,  where  k 
is  a  nonnegative  integer,  and  all  the  remaining  N  —  s  ■  k  eigenvalnes  must  be  equal  to  0. 

We  used  the  procedure  described  in  the  hrst  three  stages  of  the  algorithm  in  Section 

4  to  hnd  the  number  of  eigenvalues  equal  to  0  for  all  values  of  s,  3  <  s  <  20,  and  all 

values  of  n,  1  <  n  <  26,  where  s  is  the  number  of  sides  of  the  regular  polygon  and  n  is 
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Figure  6.  Symmetric  quadrature  of  order  32  on  the  triangle. 


the  maximum  degree  of  the  polynomials  in  the  space  U  =  7 n-  The  results  are  listed  in 
Table  5. 

We  note  that  the  entries  in  the  table  seem  to  follow  a  certain  pattern  and,  in  fact, 
they  can  be  described  by  the  following  formulae. 


Conjecture  5.1.  Let  s  be  the  number  of  sides  of  the  regular  polygon  X  centered  at  0  and 
having  all  the  vertices  on  the  unit  circle,  and  let  U  =  7n  be  the  space  of  all  polynomials 
on  X  of  total  degree  at  most  n.  We  denote  by  A(n,  s)  the  number  of  eigenvalues  equal 
to  0  of  the  operator  Pjj  oTjj,  where  the  operator  Pjj  o  Tjj  is  defined  as  in  the  statement 
of  Theorem  3.2. 

Let  us  denote  by  §p  the  sum  of  the  first  p  positive  integers: 

(5.2)  Sp  =  l  +  2  +  ...+p  =  p[p  +  l)/2  , 


with  the  convention  that  S_i  =  §o  =  0.  Then,  according  to  the  parity  of  the  number  of 
sides  s,  one  of  the  following  formulae  holds. 

If  s  =  2k  +  1,  where  k  is  a  positive  integer,  then 


(5.3) 


A(n,  s) 


§n(mods)+i,  if  n  {mod  s)  <  k  -  1 
§s-n  (mod  s)-2,  if  n  (mod  s)>k. 


If  s  =  2k,  where  k  is  an  integer  greater  than  1,  then 


(5.4)  A(n,  s) 


Sn(mods)+i  +  |(n-n(mod  s)),  ifn{mod  .s)  <k-l 

§s-n (mod  s)-2  +  ^(n  -  n  (mod  s)  +  s) ,  ifn  (mod  s)  >  k  . 


It  can  be  easily  checked  that  the  conjecture  holds  for  all  values  of  A(n,  s)  listed  in 

Table  5.  A  proof  for  the  general  case  is  beyond  the  scope  of  this  paper. 
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Constructing  quadratures  in  the  cases  when  the  operator  P  o  T  has  multiple  zero 
eigenvalues  is  more  involved  because  the  condition  number  of  the  interpolation  matrix 
(5.1)  explodes.  In  several  of  these  cases,  Newton’s  method  fails  to  hnd  quadratures 
having  positive  weights.  We  illustrate  the  situation  in  Table  6.  A  dash  denotes  that 
no  quadrature  having  positive  weights  was  found,  while  a  numbered  entry  denotes  that 
quadratures  having  positive  weights  were  found  and  the  numbered  entry  is  the  highest 
order  of  such  a  quadrature. 
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Table  5.  The  number  of  eigenvalues  equal  to  zero  for  all  the  values  of  the 
number  of  sides  from  3  to  20  (listed  in  the  top  row)  and  of  the  maximum 
degree  of  polynomials  from  1  to  26  (listed  in  the  hrst  column). 
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14 

Decagon 

- 

- 

- 
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9 

11 

13 

11 

Table  6.  Quadratures  on  regular  polygons  with  at  least  4  sides.  The 
top  row  lists  the  initial  maximum  degree  of  polynomials  n.  A  numbered 
entry  denotes  the  maximum  polynomial  degree  m  for  which  a  quadrature 
having  positive  weights  was  found.  A  dash  denotes  that  no  quadrature 
having  positive  weights  was  found. 
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